
function SimulatedMoments  = v11_Moments_all(ParametersIn);

% Parameters to be estimated
gamma         =    ParametersIn(1);
alpha         =    ParametersIn(2);
c1            =    ParametersIn(3);
c2            =    ParametersIn(4);
b1            =    ParametersIn(5);
b2            =    ParametersIn(6);

z             =    3.6;

%Set parameters
p_a = 1;
p_b = 0;
p_c = 0.5;

% True reporting rates
q_a = 0.900;
q_b = 0.409;
q_c = 0.507;
q_d = 0.234;

%Perceived reporting rates (March 24)
q_hat_a = 0.88;
q_hat_b = 0.59;
q_hat_c = 0.77;
q_hat_d = 0.23;

%Truncated true reporting costs
k_a =  9.690 ;
k_b =  4.065 ;
k_c = 5.196 ;
k_d = 2.054 ;

%Truncated perceived reporting costs
k_hat_a = 9.461 ;
k_hat_b = 6.149 ;
k_hat_c =  8.205 ;
k_hat_d =  2.007  ;


m = 20 ;



% PS effort  if s=0

num2_a =  gamma *  ((q_hat_a*m*(1-p_a))+b2);
num2_b =  gamma *  ((q_hat_b*m*(1-p_b))+b2);
num2_c =  gamma *  ((q_hat_c*m*(1-p_c))+b2);
num2_d =  gamma *  (b2);

e2n_a = num2_a / (2*c2);
e2n_b = num2_b / (2*c2);
e2n_c = num2_c / (2*c2);
e2n_d = num2_d / (2*c2);


% CHW effort if s=0


num1_a = alpha*( (q_a*((m*p_a)-k_a)) +b1);
num1_b = alpha*( (q_b*((m*p_b)-k_b)) +b1);
num1_c = alpha*( (q_c*((m*p_c)-k_c)) +b1);
num1_d = alpha*( (q_d*(-k_d))        +b1);


e1n_a = num1_a / (2*c1);
e1n_b = num1_b / (2*c1);
e1n_c = num1_c / (2*c1);
e1n_d = num1_d / (2*c1);


% PS effort if s>0

eta1_a = b2 + (b1*z)  + ( q_hat_a* ( (m* (1- (p_a*(1-z)))) - (k_hat_a*z) ) );
eta1_b = b2 + (b1*z)  + ( q_hat_b* ( (m* (1- (p_b*(1-z)))) - (k_hat_b*z) ) );
eta1_c = b2 + (b1*z)  + ( q_hat_c* ( (m* (1- (p_c*(1-z)))) - (k_hat_c*z) ) );
eta1_d = b2 + (b1*z)  - ( q_hat_d*k_hat_d*z );

eta2_a = b2 - (b1*z)  + ( q_hat_a*( (m* (1- (p_a*(1+z)))) + (k_hat_a*z) ) );
eta2_b = b2 - (b1*z)  + ( q_hat_b*( (m* (1- (p_b*(1+z)))) + (k_hat_b*z) ) );
eta2_c = b2 - (b1*z)  + ( q_hat_c*( (m* (1- (p_c*(1+z)))) + (k_hat_c*z) ) );
eta2_d = b2 - (b1*z)  + ( q_hat_d* k_hat_d*z );

dens = 2*((2*c2*(alpha^2))-(z*c1*(gamma^2)));


e2s_a = ((alpha^2)*gamma*eta1_a) / dens;
e2s_b = ((alpha^2)*gamma*eta1_b) / dens;
e2s_c = ((alpha^2)*gamma*eta1_c) / dens;
e2s_d = ((alpha^2)*gamma*eta1_d) / dens;


% CHW effort if s>0



term1_a = ( alpha * ( (q_a*( (m*p_a) -k_a)) + b1 ) )/(2*c1);
term1_b = ( alpha * ( (q_b*( (m*p_b) -k_b)) + b1 ) )/(2*c1);
term1_c = ( alpha * ( (q_c*( (m*p_c) -k_c)) + b1 ) )/(2*c1);
term1_d = ( alpha * ( (q_d*(  -k_d)) + b1 ) )/(2*c1);


obj1 = (2*(alpha^2)*c2) - ((gamma^2)*z*c1);

tempnum_a = (c2*eta2_a*(alpha^2))- ((gamma^2)*c1*z*(b2+(q_hat_a*m*(1-p_a))));
tempnum_b = (c2*eta2_b*(alpha^2))- ((gamma^2)*c1*z*(b2+(q_hat_b*m*(1-p_b))));
tempnum_c = (c2*eta2_c*(alpha^2))- ((gamma^2)*c1*z*(b2+(q_hat_c*m*(1-p_c))));
tempnum_d = (c2*eta2_d*(alpha^2))- ((gamma^2)*c1*z*(b2));

term2_a = tempnum_a /(z*q_hat_a*obj1);
term2_b = tempnum_b /(z*q_hat_b*obj1);
term2_c = tempnum_c /(z*q_hat_c*obj1);
term2_d = tempnum_d /(z*q_hat_d*obj1);


e1s_a = term1_a +(term2_a * ((alpha*q_a)/(2*c1)) );
e1s_b = term1_b +(term2_b * ((alpha*q_b)/(2*c1)) );
e1s_c = term1_c +(term2_c * ((alpha*q_c)/(2*c1)) );
e1s_d = term1_d +(term2_d * ((alpha*q_d)/(2*c1)) );




% Calculate condition for s>0 

eta3_a = c2*(alpha^2) * ( b2 -(b1*z) + ( q_hat_a*( (m + (k_hat_a*z)) )) );
eta3_b = c2*(alpha^2) * ( b2 -(b1*z) + ( q_hat_b*( (m + (k_hat_b*z)) )) );
eta3_c = c2*(alpha^2) * ( b2 -(b1*z) + ( q_hat_c*( (m + (k_hat_c*z)) )) );

eta4_a = (gamma^2)*c1*z*((q_hat_a*m)+b2) ;
eta4_b = (gamma^2)*c1*z*((q_hat_b*m)+b2) ;
eta4_c = (gamma^2)*c1*z*((q_hat_c*m)+b2) ;

eta5_a = (q_hat_a*m)*( (c2*(alpha^2)*(1+z))- ((gamma^2)*c1*z) ) ;
eta5_b = (q_hat_b*m)*( (c2*(alpha^2)*(1+z))- ((gamma^2)*c1*z) ) ;
eta5_c = (q_hat_c*m)*( (c2*(alpha^2)*(1+z))- ((gamma^2)*c1*z) ) ;

condition_a = (eta3_a - eta4_a)/ eta5_a;
condition_b = (eta3_b - eta4_b)/ eta5_b;
condition_c = (eta3_c - eta4_c)/ eta5_c;


% Calculate output
e1_d = e1n_d;
e2_d = e2n_d;
x_d = (alpha* e1_d) + (gamma * e2_d);


if p_a <= condition_a
    e1_a = e1s_a;
    e2_a = e2s_a;
    x_a = (alpha* e1_a) + (gamma *  e2_a);
else
    e1_a = e1n_a;
    e2_a = e2n_a;
    x_a = (alpha* e1_a) + (gamma  * e2_a);

end


if p_b <= condition_b
    e1_b = e1s_b;
    e2_b = e2s_b;
    x_b = (alpha* e1_b) + (gamma * e2_b);

else
    e1_b = e1n_b;
    e2_b = e2n_b;
    x_b = (alpha* e1_b ) + (gamma  * e2_b);

end

if p_c <= condition_c
    e1_c = e1s_c;
    e2_c = e2s_c;
    x_c = (alpha* e1_c) + (gamma * e2_c);

else
    e1_c = e1n_c;
    e2_c = e2n_c;
    x_c = (alpha* e1_c) + (gamma * e2_c);

end




% Side payment
st_a = term2_a;
st_b = term2_b;
st_c = term2_c;


if p_a <= condition_a
    s_a = st_a;
else
   s_a = 0;
end

if p_b <= condition_b
    s_b = st_b;
else
   s_b = 0;
end

if p_c <= condition_c
    s_c = st_c;
else
   s_c = 0;
end

s_d =0;



t_a = s_a * x_a ;
t_b = s_b * x_b ;
t_c = s_c * x_c;

de2_a = e2_a - e2_d;
de2_b = e2_b - e2_d;
de2_c = e2_c - e2_d;

dx_a = x_a - x_d;
dx_b = x_b - x_d;
dx_c = x_c - x_d;


%SimulatedMoments    = [e1_a; e1_b; e1_c; e1_d; de2_a; de2_b; de2_c; e2_d; dx_a;  dx_b; dx_c; x_d; t_a;  t_b; t_c ; e1n_a; e1n_b; e1n_c; e1n_d; e1s_a; e1s_b; e1s_c; e1s_d];
SimulatedMoments    = [e1_a; e1_b; e1_c; e1_d; e2_a; e2_b; e2_c; e2_d; x_a;  x_b; x_c; x_d; t_a;  t_b; t_c ; e1n_a; e1n_b; e1n_c; e1n_d; e1s_a; e1s_b; e1s_c; e1s_d];



return
